


archives-ouvertes 


Report on fully coupled data assimilation in simplified 
systems with implications for Earth system reanalysis 
Arthur Vidard, Rémi Pellerej, Florian Lemarié 


» To cite this version: 


Arthur Vidard, Rémi Pellerej, Florian Lemarié. Report on fully coupled data assimilation in simplified 
systems with implications for Earth system reanalysis. [Research Report] Inria Grenoble Rhóne-Alpes. 
2017. hal-01667509 


HAL Id: hal-01667509 
https://hal.inria.fr/hal-01667509 
Submitted on 19 Dec 2017 


HAL is a multi-disciplinary open access L'archive ouverte pluridisciplinaire HAL, est 
archive for the deposit and dissemination of sci- destinée au dépót et à la diffusion de documents 
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non, 
lished or not. The documents may come from émanant des établissements d'enseignement et de 
teaching and research institutions in France or recherche français ou étrangers, des laboratoires 
abroad, or from public or private research centers. publics ou privés. 





Work Package 2: Future coupling methods 


Deliverable 2.11: 


Report on fully coupled data assimilation in simplified systems with 
implications for Earth system reanalysis. 


Type: Report 
Authors: Arthur Vidard, Rémi Pellerej and Florian Lemarié 
Reviewer(s): Matthew Martin 


Delivered: 15/11/2017 


P d 


d informatics # mathematics 


Summary 


This report describes the Inria contribution to WP2 regarding a study of coupled data assimilation algo- 
rithms applied to academic models. This task was twofold: 1. Propose and study new coupled variational 
data assimilation algorithms. 2. Create a stand-alone coupled single column model (SCM) that mimics 
the ocean atmosphere behaviour and that can be used to validate the algorithms proposed in 1. 


Regarding subtask 1, a focus was made on ways to explicitly account for model coupling in the varia- 
tional optimisation problem, either as a strong contraint or as a weak constraint, or as a combination of 
both. It is then applied to both linear and non linear coupled problems, and leads to the conclusion that 
it does bring some noticeable benefit, but at a cost, both in time of development and computing time. 
The cost-benefit ratio has therefore to be studied for each given application. 


The coupled SCM has been developed in Fortran and interfaced with the OOPS data assimilation 
framework. It is documented in this report and a reference test case is given, so that it can be reused by 
partners of the project and beyond. 


1 Introduction 


In the context of operational meteorology and oceanography, forecast skills heavily rely on proper com- 
bination of model prediction and available observations via data assimilation techniques. Historically, 
numerical weather prediction is made separately for the ocean and the atmosphere in an uncoupled way. 
However, in recent years, fully coupled ocean-atmosphere models are increasingly used in operational 
centres to improve the reliability of seasonal forecasts and tropical cyclones predictions and to improve 
reanalyses. For coupled problems, the use of separated data assimilation schemes in each medium is not 
satisfactory since the result of such assimilation process is generally inconsistent across the interface, 
thus leading to unacceptable artefacts (Mulholland et al. 2015). Hence, there is a strong need for adapt- 
ing existing data assimilation techniques to the coupled framework, as presented in Smith et al. (2015). 
In that respect, ERACLIM2 is an important milestone, with the implementation and major application of 
the CERA algorithm. Task 2.11 aims at using coupled data assimilation as an opportunity to improve the 
coupling mathematical consistency of the forecast coupled system. In this report, three classes of data 
assimilation algorithms, based on variational data assimilation techniques (Le Dimet & Talagrand 1986), 
are presented and applied to single column coupled problems. Reference of fully coupled solutions are 
obtained through an iterative Schwarz domain decomposition method (Gander 2008). The aim of the 
proposed methods is to properly take into account the coupling in the assimilation process in order to 
obtain a coupled solution close to the observations while satisfying the physical conditions across the 
air-sea interface. The paper is organised as follows. The model problem and coupling strategy are given 
in Sec. 2. In Sec. 3 we briefly recall some theoretical aspects of variational data assimilation techniques, 
and we introduce and discuss three algorithms to solve coupled constrained minimisation problems. 
Comparaison with the CERA system are also presented. The performance of the proposed schemes are 
illustrated by numerical experiments in Sec. 5 in a linear case and in Sec. 6 in a non-linear single column 
ocean-atmosphere model. We decided to move technical description of our single column model and 
discussion about convergence of the data assimilation schemes in appendices, but they are an important 
part of the report nonetheless. 


2 Model problem and coupling strategy 


Numerical difficulties induced by air-sea coupling mostly come from vertical interactions, so we restrict 
our study on single column models; however most of it can be extended to 3D systems without major 
theoretical difficulties. We consider a problem defined on a bounded set O CR. Q is decomposed into 
two non-overlapping subdomains Q4 and Qo with an interface [ = {z = 0}. A model is defined on each 
space-time domain O5 x [0, T] (6 = a or o) thanks to a differential operator Z5 which acts on the variable 
ug. The problem is to couple the two models at their interface I’. To do so, we introduce the air-sea flux 
Foa and interface operators @g. Those operators must be chosen to satisfy the required consistency on 
T. 

Omitting the external boundary conditions, the equations driving the coupled column system can 
be summarised as: 


LalUa) = fa on Qa x Tw Lo(Uo) = fo on Qo x Tw 
Ua(z, 0) = uo(z) ZEQa Uo(z,0) = Uo(Z) ZEQo (1) 
Cala) = FoalUa, Uo) onT x Tw ColUo) = Foala, Uo) onT x Tw 


Where Tw = [0, T], and fp € I?(0, T; I? (Qg)) are given right-hand sides. 

In the vast majority of models, at least part of the vertical equations are solved using an implicit 
scheme, meaning that, in order to get consistency at the interface (@, (ua) = @o(Uo)), one needs to solve 
the equations on the whole column at once, which is impracticable in an air-sea context. Moreover, time 
discretisation being significantly different in the atmosphere and the ocean the meaning of this sought 
equality may not be obvious. 

Most of the time this difficulty is overcome by using so-called asynchronous coupling where the flux 
seen by the atmosphere model is computed using the (possibly averaged) ocean state from the previous 
time interval (see Figure 1), sacrificing the flux consistency on the altar of practicability. 








Figure 1: Schematics of asynchronous coupling 


Instead, Lemarié et al. (20134, D) propose to use a global-in-time Schwarz algorithm (a.k.a. Schwarz 
Waveform Relaxation, SWR see Gander (2008) for a review) to solve the corresponding coupling problem. 
This method consists in solving iteratively each model on their respective space-time subdomain using 
the interface conditions on Il computed during the previous iteration. It can be seen as iterations of 
asynchronous coupling until convergence. 

For a given initial condition uo € H! (Q4 U Qo) and first-guess u? (0, t), the corresponding coupling 
algorithm reads 


Lalu”) = fa onQax Tw [ Fo(us) = fo on Q, x Tw 
uk (z,0) = uo(z) ZEQa 4 uk(z,0) = ug(z) z € Qo (2) 
Caluk) = Foaluk ug) onPxTw | Co(uk)=Foaluk,uk) onl x Tw 


where k is the iteration number. At convergence, this algorithm provides a mathematically strongly cou- 
pled solution which satisfies 64 (Ua) = &o(uo) on T x Tw. The convergence speed of the method greatly 
depends on the choice for @ operators, and the choice of the first-guess u9(0, t). For k = 1 and if u? is 
the ocean state at previous time, then SWR is equivalent to asynchronous coupling. 

Note that this is the sequential version of the algorithm, a parallel version can be obtained using 
«6, (uE) = Z,«(uE-!, uE) instead, as interface condition for the ocean. Both systems can then be solved 
in parallel, but it is likely to require more iterations to converge. Whether the increase in parallelism 
compensates the degraded convergence depends on application and hardware configuration. 

It is difficult to advocate for SWR for operational ocean-atmosphere coupled systems, since it re- 
quires to run several instances of the model over the same time widow, significantly increasing the com- 
putational burden. However coupled data assimilation is an opportunity to improve the flux consistency 
and SWR methods can be used both as a reference and an inspiration to adapt data assimilation tech- 
niques. 


3 Dataassimilation 


Let us now suppose that some discrete estimates y of the solution to problem (1) are available over an 
irregular set of points in the interval O x Tw. In this context we are interested in using a data assimilation 
(DA) procedure to account for this additional source of information. For the present study we use the 
variational methods of DA, based on optimal control theory. Our aim is to evaluate a set of parameter xo, 
including for instance the initial condition uo of problem (1), through the minimisation ofa cost function 
J (xo) (xo is the control vector) which quantifies in some sense the misfit between the observations y and 
the model prediction. This minimisation requires the gradient of /(xo), which can be computed using 
adjoint methods (Le Dimet & Talagrand 1986). 

The cost function is generally composed of two terms: the observation term J,, which penalises the 
misfit between model trajectory and observations y} of the system, and the background term Jp which 


penalises the departure to a given reference estimate x^. 
D ui b À 0)Tp-l o 
J(x9) = (xo x ] B (xo x X (Ha (Ma, (9) y?) Rr, (Hy, (Mi; Ko) — y?) (3) 
ee i-0 
Jb ——— M 


Jo 


Where R; is the observation error covariance matrix associated with observation at time f; and B is the 
background error covariance matrix. Mr (xo) is the solution of the numerical model equations (e.g. the 
one described in equation 2) integrated from time f, to time f;, and JJ, is the observation operator that 
maps the model state space onto the observation space at time f;. 

A very important aspect for coupled data assimilation is the explicit description or air-sea error corre- 
lations (i.e. the off-diagonal blocks ofthe B matrix), but this is not the topic ofthis report and is addressed 
in other tasks of WP2 (Smith et al. 2017, for instance). The main focus here is on coupling consistency, 
so, in the following we assume that: 


Ba 0 |. : 
e B= | ^ B | is the coupled system background error covariance and B, and B, are the at- 
oO 


mosphere and ocean background error covariances matrices respectively (no explicit cross media 
error covariances are considered) 


e A()= is 0 P. 9 is the observation operator whose components are #,(.) for atmospheric 
oC 


and f (.) for oceanic observations. 


* The numerical coupled model can be written with a slight abuse of notation as 


(Mat) Mago) 
Wer dat) 0 


It includes both atmosphere (.#,(.)) and ocean (.#,(.)) as well as coupling (Maol), Moal.)) com- 
ponents 


3.1 Uncoupled variational data assimilation 


Historically the assimilation has been performed separately on the ocean and atmosphere. In a varia- 
tional data assimilation context, it means that for 6 = a and f = o, the control vector is restricted to sub- 


domain Qg and is noted xo,5 = Uolzeq,- The optimal control problem amounts to find x p the analysed 


state, which best fit observations y and a previous estimate of the initial state x) called the background. 
Noting #;,(.) the observation operator that goes from model space to observations space at time t; and 
Xg = ug the state vector, each cost function to minimise reads (for 6 = a and $ = o) 


Jg(Xgg) = (xso -x) "Bj (xso -x) «Y (Hot (Mp, 1; (%,0)) -Y5«] mo (765. (Apr 8,9)) -Y&«] (4) 


i- 
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where Rg is the covariance matrix associated to observation errors and Bg is the background error co- 
variance matrix. ./g, ;,(Xo)) is the solution of the model integration from 0 to t; starting from initial condi- 
tion x$, for B = a and B = o independently. Obviously, if the DA process is done separately on each sub- 
domain (with prescribed boundary conditions on the interface I), the initial condition up = (X7 0° X0,0) T 
obtained on Q does not satisfy the interface conditions, hence uo € H! (Q) and well-posedness of the 
coupled problem is no longer guaranteed. In practice this type of imbalance in the initial condition can 
severely damage the forecast skills of coupled models (Mulholland et al. 2015). 


3.2 Toward a coupled variational data assimilation 


Deriving data assimilation methods able to properly account for the coupling is therefore an important 
matter. This section aims at providing methods leading to a solution close to the observations while sat- 
isfying the interface conditions on T; or at least a weak form of it. The key properties of those algorithms 
are summarised in Tab. 1. 


Fully Coupled Method (FCM) 


A first possibility is to consider a monolithic view of the problem by ignoring the presence of an interface 
in the assimilation process. In this case the state vector is xg = uo(z), z € Q and for each model integration 
we iterate the models on Q, and Q, either by asynchronous coupling, till convergence of the Schwarz 
algorithm or whatever coupling method is used. 
The cost function for the FCM is 
T a qs 
Jrcu 0) = (xo - x^). B7 (xo - x^) + X (26 C460) —y5) R; (Hu (Mr, 00) y?) (5) 


i- 
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where x(t) = (ualt), uo(t)) T [tcan readily be seen that cost function (5) is identical to the cost function 
we would use for an uncoupled problem defined on Q. The solution provided by this approach is as 
coupled as the forecast system .#(.). It also allows to relax the block-diagonal nature of B and #(.), 
but for the sake of comparison with other algorithms it is not considered here. If a SWR algorithm is 
used for the computation of x(t), first-guess u? in (2) is updated after each minimisation iteration with 
the converged solution obtained during the previous model integration. The Schwarz algorithm then 
converges more rapidly over the minimisation iteration. Note that the FCM requires the adjoint of the 
strongly coupled model (2) which can be tedious to derive. The main drawback of this method with SWR 
coupling is that it possibly requires a very large number of Schwarz iterations since it systematically 
iterates till convergence. 


Partially Coupled Method (PCM) 


In this variant only one part of the coupling process is accounted for when computing x(t) in the cost 
function. It is aimed at improving the coupling consistency of the solution through data assimilation by 
adding a penalty term in the cost function. For instance one can propose to truncate the Schwarz itera- 
tions in the direct and adjoint model after kmax iterations, with kmax < kevg. Because we do not iterate till 
convergence, the coupled solution does not strictly satisfy the interface consistency. The equivalent ap- 
proach in asynchronous coupling (kmax = 1) would be to use one-way coupling to compute the trajectory 
and to aim at 'promoting' it to two way through a penalty term in the cost function. 

As proposed by Gejadze & Monnier (2007) in the context of river hydraulics, a convenient way to 
propagate the information from one subdomain to the other during the minimisation iterations is to 
use an extended cost function which includes the misfit in the interface conditions. The idea behind 
this approach is to enforce a weak coupling within the minimisation iterations. The control vector xo = 
(uo(z), u? (0, 2)? now includes the first- guess on the interface and the cost function reads 


N 
Ipc (Xo) = J^ (xo) + X (Ha C41 ""*(xo)) — y) Rz! (265 CAT UNE (xo) — 2) + 5 Ko) (6) 
i=0 


where 4/7" (x0)) = (uv (1), uv (1:))T and 


J Go) = liés lUa (0, D) — (90, Do, T) (7) 


Unlike FCM, the first-guess for the interface is part of the control vector here, but this method still 
requires part of the adjoint of the coupling. Note that since the first-guess u? is updated at the end of 
each minimisation iteration, we can expect that we will converge toward a good approximation of the 


SWR solution. 


Weakly Coupled Method (WCM) 


The last possibility we propose to investigate is to suppress the coupling operators in the model and rely 
solely on the minimisation iterations to weakly couple the two models. This approach only requires the 
adjoint of each individual model but not the adjoint of the coupling as for the previous algorithms. The 
control vector is Xo = (X4,0, Xo,0) T with Xo,p = (Uol zea ges ub (0, £)). The corresponding cost function is 


Twom (Xo) = JP(x4,9) + JP (6,0) + J2(x4,0) + J2 (5,9) + J (Xo) 


It is straightforward to see that this algorithm provides only a weakly coupled solution. One model in- 
tegration is performed (possibly in parallel for ocean and atmosphere) with boundary conditions on T 
provided by the term us (0, t) taken from the control vector. 


Table 1 summarises different aspects of the three proposed algorithms. In the 'adjoint of the cou- 
pling’ column, online means that it is required along with the adjoint model during the VJ’ computa- 
tion, while offline means that it is only required during V/? computation, so that contributions to VJ’ 
from ocean and atmosphere can be computed in parallel. 








Algo Control # of extended Adjoint of the coupling Coupling 
vector coupling cost 
iterations function 
FCM (uo(z)) Kcvg no online strong 
PCM  (ug(z) u9)* Kmax yes partially online/offline ~strong 
WCM  (uo(2,u9, u2)? 0 yes offline weak 





Table 1: Overview of the properties of the coupled variational DA methods described in Sec. 3.2. Nota- 
tions are consistent with those introduced in the text. 


4 Incremental formulation and link to the CERA system 


In classical 4D-Var, due to non linearities in 4 and #, minimising efficiently (4) is not straightfor- 
ward. Common practice is to use the so called incremental 4D-Var approach (a.k.a. Gauss-Newton in 
the optimisation community) where the original problem is solved through successive minimisations of 
quadratic cost functions (inner loops) 





Ex a l-1 N 

J'(6x') = [ox +Y sx) B? [ox +Y os +) (Hi ME ‘6x! - a; Re (ni ME ax -aL! e) 
k=1 k=1 i=0 

with d? = yr, — Hn (Mn (Ko)) and d} =y}, -Hr (Mr Xo+ a óx!)) 

and M} (resp. H} ) being the tangent linear operator of Mr, (resp #%;,) differentiated around Xo + 
vis 6x! 

Non linearities are therefore accounted for through the re-linearisation of the M' and H! operators 
and the computation of the innovation vectors d/. Under some regularity hypotheses, such algorithm is 
known to converge toward the solution of the original problem (see appendix C for more details). 

FCM being a direct transposition of classical 4D-Var, its incremental variant uses the same inner loop 
formulation as equation 8. For PCM and WCM, both J} and J : can be obtained similarly as equation 8 
as well, and the inner loop coupling penalty term reads 


J! = yllCaðu} (0, [) — C,6ul(0, D) satu, (0, t)) — «e, (ul! (0, t)) loa 


C, and C, being the tangent linear operators of 64 and €, respectively. 
The CERA system uses a different formulation compared to regular incremental 4Dvar, indeed it aims 
at solving the FCM problem by successive inner loop approximations that are uncoupled, i.e. CERA uses 


l 
MI) = pu E in the outer iterations, but M! = | Ps d 1 
speaking, CERA can be seen as a parallel Schwarz algorithm for solving FCM minimisation, with the 
important limitation that it cannot cope with air-sea correlation in B. Likewise, PCM with no J; term 
(y = 0) can be seen as a sequential Schwarz algorithm for solving FCM. 


Convergence properties of algorithms presented in this report are discussed in appendix C. 


| in the inner loops. Roughly 


5 Application to a 1D linear diffusion problem 


In this section, previous algorithms are applied on a 1D diffusion problem. We, thus, consider for both 
part of the system Zp = 0, vgO? in (2), with B = a and B = o and with v4 Z v, the diffusion coefficients in 


each subdomain. The computational domain is Q =] — La, Lol with Lg, Lo € R**. We choose the interface 
. - a ; 0, 0 1 0 
operators on T to obtain Dirichlet-Neumann condition, i.e. €, = | e = 1 | and 6, = | 0 vô | . We 
0 Z 


consider the analytical solution Up» and the corresponding right hand side fg = u^, of the coupled 
problem on each subdomain as : 


us(z,t) = —e “b 434 cos^| — on Og x Tw. (9) 
p 4 T 


where Up = 20 °C and 1 = 22 h. Note £4v, = EoVa is required to ensure the proper regularity of the 
coupled solution across the interface I'. To satisfy this constraint we choose £4 = 4 km, €, = 0.4 km, 
va = 1 m?/s, v; = 0.1 m?/s. The model problem (2) is discretized using a backward Euler scheme in time 
and a second-order scheme in space. The resolution in each subdomain is Az = 20 m with Lg = Lo = 1 km 
and the time-step is At = 180 s. The total simulation time is T = 12 h, which is also the size of the SWR 
window (Tw in equation 2). The latter implies that at least 2 iterations of SWR are necessary to get some 
coupling in the models integration. 

For the assimilation experiments, we consider true-state x! to be the analytical solution while back- 
ground x? corresponds to the solution obtained with a biased initial condition. In general, the Schwarz 
algorithm converges in keyg = 50 iterations with a tolerance € = 1078. Some observations y of the true- 
state are generated such that y = H(x’), with H the observation operator. The observation and back- 
ground errors covariance matrices are considered diagonal such that R = 10 Id and B = 100 Id. For the 
extended cost function we consider different values of y. All the minimisation are done until conver- 
gence of a conjugate gradient algorithm with a stopping criterion || V (xo) lloo« 107°. 


Single column observation experiment 


For our experiments, we consider that observations are available in © \ {I} only at the end of the time- 
window (i.e. at t = T). In this case, the results obtained for different assimilation schemes are reported 
in table 2 where the performance of each scheme is presented both in terms of number of minimisation 
and computational cost. The latter being given relative to that of uncoupled data assimilation. Possible 
parallel aspects are not accounted for in thes measure of computational cost. To evaluate the strength 
of the coupling we define an interface imbalance indicator which corresponds to the value of J; at the 
end of the DA process, with y = 1. Values of J* close to zero indicate a strongly coupled analysed state. 


In table 2, a root mean square error (RMSE) defined as 4/ E ((x^ — x/)?) on Q x Tw is also used to evaluate 
how much the analysed state is close to the true-state. 





Experiment names start with the algorithm name, then, if relevant, the number of Schwarz iteration 
(F stands for full convergence) and finally, if relevant, Jsy value. From table 2, we can first note that the 




















Algo Y  Kmax #ofminimi- Computing Interface RMSE 
sation cost imbalance in°C 

iterations (relative to indicator 

uncoupled) 

FCM-F - kcvg 14 28.5 3.10 1? 0.221 
FCM-2 - 2 14 1.3 9.866 0.415 
PCM-F TEE kog 14 28.5 3.102 0.221 
PCM-5 0 5 14 3.1 4.107? 0.220 
PCM-2 0 2 14 1.3 9.87 0.415 
PCM-2-Js0.1 0.1 2 184 15.8 2.1074 0.216 
PCM-1-Js0.1 0.1 1 117 5 6. 107? 0.217 
WCM-Js1. 1.0 0 365 16.3 1.1078 0.286 
WOCM -Js0.1 0.1 0 396 16.9 1.106 0.228 
WCM-Js0.01 0.01 0 390 16.6 1.1074 0.226 
Uncoupled 0 0 22 1 13700 8.338 





Table 2: Results obtained for the three coupled variational DA methods described in Sec. 3.2 with several 
settings. Observations are available in Q \ {I} at the end of the time-window. 


FCM algorithm requires few minimisation iterations to obtain a low RMSE value and a strongly coupled 
analysed state (J? ~ 10 1?). A drawback of this approach is a high computational cost (almost 30 times 
that of uncoupled). PCM-F (resp -2) only differs from FCM-F (resp -2) by the first guess of the inter- 
face within the control vector, and gives pretty much the same results as FCM-F (resp -2), for the same 
computing cost. This shows that without /* adding the interface in the control has little effect. 

Since in the other PCM approaches, coupling iterations are truncated and first-guess u? is part of 
the control vector, we expect a reduced computational cost compared to FCM-E If y = 0 (no J; term, 
PCM-2 and 5) the same number of iterations are required for the minimisation to converge reducing 
dramatically the cost over FCM, but at the expense of a lower quality analysis. On the one hand de- 
creasing the value of kmax increases the number of minimisation iterations. Indeed, going to Schwarz 
convergence (kmax = keyg) procures the best model solution, it then needs few minimisation iterations. 
However, for the next iteration, the background interface is given by the control vector rather than the 
previous converged estimate; therefore it requires again numerous Schwarz iterations. On the other 
hand, by reducing the kmax value, the number of Schwarz iterations is reduced and the update of first- 
guess more significant. However the coupling quality is affected and this leads to a slower minimisation 
convergence. Here, a good compromise is to choose kmax = 5. 

When J; term is activated (PCM-2-Js0.1 and -1-Js0.1) the number of iterations required for the min- 
imisation to converge rises significantly as it is usually the case when one adds constraint to the cost 
function. However, computing cost is still significantly lower than that of FCM-F while reaching similar 
quality analysis. Smaller values of kmax provide a faster convergence of the algorithm. With kmax = 1, 
which corresponds to a one-way coupling, it requires only 5 times the cost of the uncoupled DA to pro- 
vide a good approximate of the strongly coupled solution (J = 6.87 10 ?, RMSE = 0.217°C). For kmax > 1 
two ways interactions are accounted for in the model integration, making interactions between the con- 


trol vector components more intricate and therefore damaging the convergence properties. 

By considering uncoupled models in the WCM algorithm, the convergence property of the minimi- 
sation is severly damaged, as can be seen in figure 2 This is true whatever the choice of y to balance /? 
and J? in the cost function and leads to a significant increase of computing cost compared to PCM. The 
analysed state shows a larger interface imbalance indicator compared to FCM and the best PCMs, which 
confirms that WCM provides a weakly coupled solution, but is significantly better than the uncoupled 
DA in that respect. The RMSE level is similar to that of FCM and most PCMs though. 
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Figure 2: J° evolution for the different schemes, along minimisation iterations (left panel) and corre- 
sponding computing cost (right panel, unit: uncoupled minimisation cost) 


Table 2 presents results at convergence of the minimisation. In realistic ocean-atmosphere coupling 
such convergence is generally not affordable so rate of convergence may be more important than the 
actual minimum. Table 3 shows the same results as 2 for the best PCM and WCM limiting their cost to 
that of uncoupled minimisation. 














Algo Y  Kmax #ofminimi- Computing Interface RMSE 
sation cost imbalance in °C 
iterations (relative to indicator 
uncoupled) 
PCM-1-Js0.1 0.1 1 22 1 5.10 4 0.353 
WCM -Js0.1 0.1 0 22 1 2980 2.51 
Uncoupled 0 1 22 1 13700 8.338 





Table 3: Same as previous table but limiting the number of minimisation iterations to 24 


Both approaches outperform the uncoupled DA by far, PCM-1-Js0.1 being even quite close to its 
converged state, with less than a third of its original cost. For WCM to be a viable option, in that con- 
text, some substantial work has to be performed on preconditioning, in order to significantly improve 
convergence. 

In the linear case, a CERA equivalent can be defined using outer and inner loops, with only the outer 
loop being coupled. Here the incremental approach is used to account for the coupling instead of the 


10 


non linearities. Any strength of coupling can be used in the outer iterations. Results are given in ta- 
ble 4 for fully converged coupling and 2 iterations of Schwarz (corresponding FCM results are recalled). 
Both CERA-F and CERA-2 give similar results to their FCM counterparts with a much higher iteration 
count, but a smaller computational cost. In CERA-2 the coupling Schwarz iterative scheme restart from 
the same initial interface condition at each outer iteration. CERA-2-save on the other hand, reuses the 
output of the Schwartz algorithm from the previous outer iteration, which allows to really improve the 
strength of the coupling for very little extra-cost. This improvement is really related to the way the cou- 
pling is done in this case, though. Indeed since the SWR window is the same as the assimilation window, 
the background of the interface has a strong impact when using low SWR number of iteration. 











Algo Y kmax #ofminimi- Computing Interface RMSE 
sation cost imbalance in °C 
iterations (relative to indicator 
uncoupled) 

FCM-F - kwg 14 28.5 3107 0.221 
FCM-2 - 2 14 1.3 9.866 0.415 
CERA-F  - kog 188 15.5 4.10 7? 0.271 

CERA-2 = 2 190 8.95 10.37 0.49 
CERA-2-save — 2 194 9.13 8.107? 0.270 





Table 4: Same as previous table but for FCM and CERA experiments 


Fig3 presents the evolution of outer and inner cost functions during minimisation for CERA-E CERA- 
2 and CERA-2-save. It shows that CERA-2-save benefits from previous outer iteration to improve its cou- 
pling convergence and even outperform FCM-2. On the downside, figure 3 right panel shows erratic 
inner minimisation behaviour with some oscillation around the optimum making the choice of stop- 
ping criterion complicated. It also shows a degradation of the convergence properties of CERA-2-save 
successive inner loops (growing offset between vertical dashed lines). 

It is difficult to draw a definitive conclusion from such a simple test case. The next section goes 
a step further toward realistic applications and present preliminary results on a more complex ocean- 
atmosphere single column model. 
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Figure 3: left panel: outer cost function values for different flavours of CERA (in red) and corresponding 
FCM optimum value (horizontal green lines). Right panel: inner cost function values for CERA-F and 
CERA-2-save vertical dashed lines mark outer iterations. 


6 Application to a single column ocean-atmospheric boundary layer model 


In order to perform further tests of our algorith- z 
mic developments, we set up a more realistic single ha 
column model of a column of the ocean coupled 
with a boundary layer of the atmosphere. The at- Atmosphere 
mosphere is defined, as presented opposite, on do- Qa 
main Q, = [0, ha] and the ocean on Qo = [-ho,0]. 

The surface boundary layer is located between z* CLSA 
in the atmosphere and z in the ocean. Ocean at- Interface r 0 E CLSO 
mosphere interface being represented by r. The Z 

time domain is [0, T] with T > 0s. We can then 

express evolution equations for velocities u and Ocean 

tracers t (temperature and humidity for the atmo- Oo 

sphere, temperature and salinity for the ocean) of -ho 

our 1D single column model : 




















Oug(z, t) ô dug(Z, t) 
TT =~ flx ugla D+ E (KRO EE) + rss 0 sur Og x [0, T] 
ôtg (z, t) Ô p Otg (z, t) 
Ty = az K; D + Ft, (z, t) sur Og x [0, T] 
ô ô 
aK (zt) 2] = pK (2) | = £7 (us, uo tato) sur T x [0, T] 
Oz T Oz F 
ar + Ot Of„— Ot; ors 
PaK; (Z) — | -poK;(z ) =—| = Foala Uo, ta, to) sur T x [0, T] 
Oz T Oz T 

















where B = a, o refer to atmosphere and ocean variables respectively. Both models use the same structure 
and differ from their forcing terms F,, their interface conditions and the computation of their turbulent 
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Figure 4: Truth (plain line) and background (dashed lines) initial conditions for atmospheric (top) and 
oceanic (bottom) quantities 


viscosity and diffusivity coefficients KÊ, and KÊ . The construction of this model from the 3D equations, 
along with the definition of air-sea fluxes 7}, are described extensively in appendix A 


6.1 Model and data assimilation settings 


The model settings used in this study are described in appendix B. The assimilation window is chosen to 
be 12h, starting from day 2, so that it includes a switch from stable to unstable regime. The SWR window 
is set to 1h, so that one-iteration only will yield a classical asynchronous coupling with exchange of fluxes 
every hour. Observation are taken every 3h and all the vertical column is observed and perturbed with a 
white noise, corresponding to specified observation errors variances in R, namely 


Tu, = 0.4 ry, = 0.0075 
ry, = 0.4 ry, = 0.0075 
“o (10) 
rg, = 0.35 rg, — 0.8 
Yg, = 0.0002 rs, — 0.055 


The background is obtained by perturbing the true state 3 days prior to the start of the assimilation 
window and propagated forward in time so that it is physically balanced. Both truth and background 
initial quantities are shown in figure 4. 


The background error correlation matrix is considered block diagonal (no cross variable covariances), 
each block representing Gaussian covariances. 
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" 0 B, 0 Ol _| 0 B, O 0 
^" |o 0 Bo, Of * |O 0 Beg O 
0 0 0 B 0 0 0 Bs, 
5 da an 
Vom, 0 of [m 0 0 
Bx,r = Ka Bx =| 0 B. 0 
a, 0 0 Bo, 0 |” ES ON 
t 
0 0 0 Bay = 
Let b; ; being the B coefficient at line i and column j. It is set to: 
k=min(|j-il, N — j+i) where N is the number of vertical levels 
tan) 
b; ; -o^exp|- (12) 
5j p | 212 
bi; = bij 
In this experiment, we also set: 
gy 5 Lu, = 50 04,5002. 15 
d s Ly, = 50 Ov, =0.02 Ly, =15 
O9, =1 Lg, = 40 O0, =1 Lg, = 50 
Oq,= 0.0007 Lg, = 20 05,702 Ls =30 
(13) 
Ouar =1 Luar =9 Ossu = 0.1  Lssu =O 
Del Lor =0 UO Lue 
Tg, = 10 Ig 50 Osst=15 Lest =0 


Ogar =0.001 Lg, =0 
Finally, the coupling penalty term is defined for PCM and WCM respectively as 


Joem D = a (Foa xb, xk, B) — Foa xb xk, 2) W (Foa xE, R) - Font ism) (19 
and 
Tivem (Xo) = a (Foa(xD, xo, À) — Foa Kax, À) TW! (Pak, Xo, 22) — Foa XaX?, À) (15) 
with a = 0.001 and W such that: 


W, 0 0 0 
0 W, 0 0 


W = 16 
0 0 Wo, 0 
0 0 0 Wr 
where 7, are wind stresses, Qnet is the net heat flux and F is the fresh water flux: 
wz; = 0.001 
wz; = 0.001 
(17) 
WQret = 1 
wr =310° 
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In order to mimic realistic settings, only two outer iterations will be performed. The inner minimisa- 
tion will be allowed to go to convergence since they are quite efficient with these settings. 


6.2 Numerical results 


Table 5 presents a similar summary as of previous section, with the notable difference that the analysis 
error is represented as improvement over background, meaning the higher the better (background being 
at 0% and truth at 100%). It is computed as the mean over physical quantities of 


lel? — e^ ||? 
lle? |? d 
where &" and e^ are background and analysis error respectively. 
In this more demanding setup, CERA does not yield the same result as its FCM counterpart. In- 
creasing the number of outer iterations (not shown) actually degrades the CERA results, hinting that the 
convergence requirement are not met and that the model used in the inner loop is too different from the 
outer one. This is aggravated in CERA-F where the outer model goes to SWR-convergence. This lack of 
outer convergence precludes CERA from fully benefiting from outer coupling 
On the other hand, accounting for the coupling during the assimilation process (through PCM or 
WCM) allows to get a reasonably good analysis, both in term of coupling and RMSE. 

















Algorithm y kmax #ofminimi- Computing Interface RMSE 
sation cost (relative imbalance improvement 

iterations to CERA) indicator (in %) 
FCM-F  - kw 26 3.8 210 74 
FCM-1 = 1 28 1.06 5 65 
CERA-F  - kog 24 14 5.81071? 24 
CERA-1 = 1 26 1 1.6 40 
PCM-1 01 1 25 0.96 4.10? 60 
WCM 0l 0 31 1.2 6. 1078 57 





Table 5: Result summary for the SCM system 


Differences are largely located at the limit of the boundary layers that CERA and, to a lesser extent, 
WCM tend to misplace. Figure 5 shows the difference in the analysed state between FCM-F and CERA-1, 
FCM-F and PCM, and FCM-F and WCM for atmospheric temperature and v-velocities along the assimi- 
lation window. For temperatures, even if CERA benefits from its less constrained optimisation problem 
to get a better temperature at the beginning of the time window, not accounting enough for the coupling 
processes quickly degrades the analysis. PCM and WCM manage a better estimation of the ABL height 
and are close to the FCM one throughout the assimilation window. One can notice that as for CERA, 
PCM benefit from more degrees of freedom and improve the analysis compared to FCM at the beginning 
of the window, and thanks to its stronger coupling, it manages to retain a good analysis. 

Differences are less striking on ocean quantities (not shown) and actually CERA-1, with a 40% im- 
provement over background is doing a reasonable job. Longer term forecast can also be affected by 
differences in the initial conditions. Figure 6 shows the evolution of sea surface velocities for four days, 
starting from the end of the assimilation window, using the same model coupling (asynchronous). All 


15 





2000 2000 2000 








1500 1500 1500 


zy (m) 


1000 0.00 


Mme p —0.08 
"T 
500 = 


1000 0.00 1000 0.00 


500 


500 pet 


“Error difference FCM-CERA (in %) 
(m) 

2 Error difference FCM-PCM (in 96) 
r (m) 

` Error difference FCM-WCM (in %) 























2000 2000 








2000 





1500 1500 1500 


1000 1000 0.00 


i 


y m 
Error difference FCM-CERA (in %) 
x (m) 


500 


500 


. Error difference FCM-PCM (in %) 
(m) 


Error difference FCM-WOM (in 96) 


-032 





























t (h) t (h) t (h) 


Figure 5: Differences in the analysed state between FCM-F and CERA-1 (left panels), FCM-F and PCM 
(middle panels), and FCM-F and WCM (right panels) for atmospheric temperature (top) and v-velocities 
(bottom) along the assimilation window 


three systems improve significantly from the background and get the right variability even though fore- 
cast from CERA analysis has a tendency to overshoot. 
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Figure 6: Forecast of SSU and SSV from FCM-E CERA-1 and WCM analysis. Dashed and plain black lines 
are background and truth evolutions respectively 


With a 20% increase in computing cost WCM brings a noticeable improvement over CERA, without 
most of the complexity of an FCM or PCM scheme. It therefore and may be a good candidate for an 
intermediary step toward fully coupled data assimilation. However this optimistic conclusion needs to 
be moderated. First Note that in this study, inner problems in CERA are solved as a single optimisation 
problem, unlike realistic application where ocean and atmosphere are solved separately. Although it 
does not change the result at convergence, doing so may affect the convergence speed of the inner prob- 
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lems and therefore increase artificially CERAs computing cost. Also experience showed that WCM was 
more difficult to properly tune and required much more trial and error to get it working efficiently than 
the other algorithms. On a simplified model it is doable, in order to go toward more complex problems, 
one will need to devise an efficient tuning strategy. 


7 Software developments 


Both linear and non linear coupled models used in this report have been developed and interfaced within 
the OOPS environment. The coupled diffusion model is made of 9000 lines of C++ including hand- 
written tangent and adjoint models and interfaces while the coupled SCM is made of 70 000 lines of For- 
tran/C++ including tangent and adjoint and interfaces. SCM tangent and adjoint models were generated 
using the automatic differentiation tool TAPENADE!. It will then be made available to the community 
and be proposed as a reference test case to be included in OOPS. 


8 Conclusion and perspectives 


Ocean-Atmosphere coupling is a complicated matter and is still a somewhat open question. In this re- 
port we mostly mentioned the time inconsistency in the fluxes exchanged by ocean and atmosphere, 
one can also add the highly parameterised nature of such exchanges and the uncertainties that are asso- 
ciated. In both cases coupled data assimilation is a challenge but also an opportunity. First by improving 
the coupling consistency of the analysis, but also, going beyond the point of this report, by providing ma- 
terials for improving coupled modelling. The latter can be achieve by studying the coupled increment 
statistics, for instance, provided that the coupling is accounted for in the data assimilation scheme. 

In this task we compare several possibilities to improve the coupling quality through variational data 
assimilation. Three algorithms are presented and their convergence properties discussed in appendix C. 
They are applied first to a simple linear diffusion coupled problem and then to a more realistic coupled 
single column model. The coupled SCM is thoroughly described in appendices A and B and is itself a 
result of the task, since it can be used as a testbed for future research in coupled data assimilation. 

One can draw a few conclusions from this preliminary study. First, CERA was probably a reason- 
able choice, with a good trade-off between low complexity and efficiency. Indeed, if both ocean and 
atmosphere data assimilation systems are available, its implementation is relatively easy (as operational 
implementation permits) and it is able, thanks to its outer loops, to account for part of the coupling. 
However, this is limited by its lack of global convergence. The latter being also observed on CERA ap- 
plied to operational OA system, where nothing is gained beyond a couple of outer iterations (but much 
is gained by the second outer iteration). Second, accounting explicitly for the coupling, either as a strong 
or a weak constraint of the variational data assimilation problem, can be of benefit. Indeed, it improves 
the overall quality of the analysis and the global (outer loop) convergence of the system. As a side prod- 
uct it gives a feedback on the coupling processes, that can be highly valuable to improve the system, as 
mentioned above. As a downside, it increases the complexity of the system, hence its development time. 
Moreover it potentially damages the local (inner loop) convergence speed, hence its computing time. 

The real conclusion of this report is that ocean and atmosphere coupled modelling and coupled data 
assimilation are both a complicated matter, and one should look into both problems at the same time in 
a consistent manner. 





lhttp://www-sop.inria.fr/tropics/tapenade.html 
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A Single column coupled model description 


This appendix describes a single column coupled ocean-atmosphere model. It is a quite valuable tool 
since it mimics the complexity of model used in a forecasting center without the tremendous compu- 
tational burden. Single column models are also widely used to study sub-mesh parameterisations in 
the ocean-atmosphere surface boundary layer without the cumbersome 3D model. These sub-mesh 
parameterisations introduce non-differentiabilities in addition to strong non-linearities. These non- 
differentiability and non-linearities can cause several problems in a process of data assimilation and the 
development of an adjoint model (Janisková, Veersé, Thépaut, Desroziers & Pouponneau 1999, Janisková, 
Thépaut & Geleyn 1999), which is why it will be very interesting to use such a model to test the robustness 
of our algorithms introduced in section 3.2. 


A.1 Modeldescription 
A.1.1 From3D primitive equations to single column 


Equations of our simple column model derive from 3D primitive equations used in atmosphere or ocean 
operational forecasting models. These are the Navier-Stokes equations describing the motion of a New- 
tonian fluid simplified by various hypotheses themselves specific to the case of atmosphere or ocean: 


* Fluid is in hydrostatic equilibrium, which consists in neglecting the vertical acceleration of the 
fluid in equations of vertical motion. 


* One neglects the vertical component of Coriolis acceleration. 


* The thin-layer hypothesis is carried out, that is to say that the fluid is contained in a very thin layer 
relative to the radius of the sphere. In practice this amounts to neglect altitude with respect to the 
earth radius. 


These equations are thoroughly studied in the literature (see McWilliams (2006) for example) 


A.l.l.a Single-column hypothesis and Reynolds decomposition 


In order to study the behaviour of the models at the air-sea interface with sub-mesh parameterisations, 
it is possible to make some additional simplifying hypotheses. These hypotheses allow us to reduce to a 
simpler case study, in 1D, while keeping the complexity of modelling parametrisation at the interface. We 
will therefore have a somewhat representative approximation of a realistic model based on primitive 3D 
equations with the simplicity and low numerical cost of a ID model. Such a model is obviously not meant 
for weather forecasting, but still behaves in a similar way than that of a 3D realistic model. Therefore, we 
can hope that the results obtained will be transposable to 3D models. 
The first hypothesis is to consider horizontal homogeneity: 


ô __ à _ 
. 3x Bye 


From there we can simplify the 3D primitive equations, leading to: 
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ou owu 
(Momentum equation) Sr =-fkxu+ ETA + v, V^u tu 
Z 
: : Ow 
(Conservation equation) E 0 (19) 


(Tracers equation) OP 
Where u = (u, v)" and w are zonal, meridional and vertical wind or sea current velocity components, f 
is the Coriolis factor and , vm and v; are molecular viscosity and diffusion coefficients respectively. 9^ 
represents the various tracers considered in the atmosphere and the ocean (air humidity q4, air and sea 
potential temperature 05 (6 = a,0) and sea salinity So). Finally Fu,y terms represent the sources and 
relaxation terms. 

Applying the Reynolds decomposition to the above equations, one gets: 








ou) | 0(w)(u) O(w'w) à » 

ax ES az t oz + vg V^ U) + Fu 

OT) OWT) OWT’) AEE 20 
age gm ge HEN (T)+ Fr ( ) 
ðw) | 

Oz E 


A second hypothesis is to neglect the vertical advection, i.e. (w) « 0 allowing for a further simpli- 
fication. The unresolved sub-grid variations (w’u’) and (w'") are represented through the turbulent 
closure scheme: 


( lu) = K. ô (u) 
Oz 
(21) 
ies 
* óz 


where Kn is the turbulent viscosity coefficient and K; the turbulent diffusivity coefficient. 


A.1.1.b Single column model equations 


A 1D single column model can then be derived i 
from the above mentioned hypotheses. First by ha 
defining the vertical domain of our model (figure 
opposite). The atmosphere is defined on domain Atmosphere 
Qa = [0, ha] and the ocean on Qo = [-h,,0]. The Qa 
surface boundary layer is located between z* in ü 
the atmosphere and z in the ocean. Ocean atmo- a CLSA 
sphere interface being represented by I'. The time InterfaceT 0 E CLSO 
domain is [0, T] with T > 0s. We can then express Z 
the equations of our 1D single column model: 
Ocean 
Qo 
—ho 
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Oug(z, t) Ó Oug(z, t) 

— T = — flex up(z, 1) + — KE = | + Fag(z,t) on Og x [0, T] 

04 p(z, t) ô f OT p(z, t) ee) 
a: = 3z (x Z az EZE t) on Qg x [0, T] 











with $ = a or o designating atmospheric or oceanic variables. Thus, both atmosphere and ocean model 
have the same structure and differences will come from RHS choices, interface conditions and compu- 
tation of viscosity and turbulent diffusivity coefficients. 


Molecular viscosity and diffusivity terms vanished because they are negligible compared to turbulent 
terms within the surface boundary layer. Outside the surface layer, turbulent coefficients can be consid- 
ered as being constant and equal to molecular coefficients and representing the free atmosphere or the 
deep ocean properly. Note also that only resolved terms (u,,;) and (2,9) are present in the equation, so 
(.) are omitted in 22 for the sake of clarity. 


In our model, we use the following forcings Fy, ,(Z, t) and F. Tao (2 É): 


V(z t) € Qao x [0, T] : 
Fu, (Z, t) = — fk x uaa, t) relaxation toward geostrophic winds u,G(z, t) 
Fu, (Z, t) = 0 


Fg, (z, t) = As(z) (04 (z, t) - Ors(z, t)) large scale relaxation toward 07 s(z, t) 





(23) 
rt : 
Fa, (2, f) = 0z Qs L penetrating solar fluxes Qs(z, t) 
PoCp 
Fq,(Z, t) = As(z) (da(z,  — qrs(z, 0) large scale relaxation toward qzs(z, t) 
s (z,t) -0 


with 1/As(z) a altitude dependent relaxation time allowing to define an atmospheric large scale nudging 
term and C, water thermal capacity. Qs(z, t) is the solar flux penetrating into the ocean, computed from 
the surface solar flux Qs, = Qs(0, t). Underscore zs stands for Large Scales, these quantities being defined 


externally. 


The model also require conditions at z = hmax and z = — Pmax, external boundaries of the domain, 
denoted 00% : 





dug t 
———(t) -0 on aag 
Oz 24) 
035 ( 
_ ext 
dE (t) = F3; on 008 
where tracers fluxes F T, are chosen as: 
006 (z,0 
Fg, =0 Fg, = ve Em 
z=-ho 
0S, (z,0) es 
Fy, =0 Fs, = v. ———— 
Oz z=—ho 
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A.2 Interface boundary conditions and air-sea coupling 


These conditions at ocean-atmosphere interface must make it possible to verify the continuity of mo- 
mentum and tracers fluxes. Both models are therefore forced at interface I' by different fluxes calculated 


from atmospheric and oceanic components, which allows for a coupling between the two media. 


A.2.1 Momentum flux continuity 


Recall that atmosphere and ocean velocities follow: 





du à 
a = — fk x (ua — Uag) + Oz (K5,0zu,) 


OU, 
Ot 


Interface conditions on ug and u, are given by: 





ô 
=-fkxu,+ E (K Oz) 


ôu ôu 
Pak gz, = Poking, 


where r is the momentum flux, defined by the bulk formula: 


=T onlx{0,7] 


T = PaCa | ua(z*) - uo(z ) || (ua(z*) -uo(z )) 


with Cg the bulk drag coefficient. 


A.2.2 Heat fluxes 


In the ocean, potential temperature 0, in (22) and (23) is defined by: 


de 2 (ue, Qs. 
dt dz\ * Oz poCy 





At ocean surface, flux continuity imposes: 


0060 | Qs, = Qnet 
` AZ PoC,  poCp 








on I x [0, T] 


where Qnet is the net heat flux, defined as: 
Qnet = Qs, + Qr, + Qrt + QE + Qu 
with: 


* Qs, surface solar flux (prescribed) 


(26) 


(27) 


(28) 


(29) 


(30) 


(31) 


* Qr (t) long wave radiative input normally coming from the atmosphere model, but prescribed 


here. 


e Qr; long wave radiative output. Considering the ocean as a black body, Qz; is given by the Stephen- 


Boltzmann law: Qr; = —ea0,(z-)^ with e = 1 ando = 5.67- 10 ?W.m ?.K 4. 


e Qg and Qz sensible and latent heat fluxes (defined later) 
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Regarding the atmosphere, temperature satisfies: 











00, O 00, 
- K? À — 32 
ôt sl Ss m s (da Ors) (32) 
At the interface: a0 Q 
| und eR Li Tx [0,T 33 
Oz pals 50 (0, 2 Gy 


with C the thermal capacity and Qy the sensible heat flux defined by: 


Qu = paC5 Cn || walZ*)) - uo(z )) Il (Aa(Z*)) — 5C )) (34) 


with Cy the bulk sensible heat transfer coefficient. 
Note that we made an additional simplification: the atmosphere is supposed to be transparent, i.e it 
does not absorb radiative fluxes, so they do not directly impact air temperature. 


A.2.3 Evaporation 


Air humidity qa is defined by 





04a 0 aa 
E K 
ôt 5l 5 az 


| * As(da — ars) (35) 


The interface condition is: 





ðqa -E 
Ke- oon rx (0,71 (36) 
ÔZ Pa 


where E is the evaporation flux defined by bulk formula: 


E = paCr || ua(z^) - uo(z ) | (ga(z*)— qo(z )) (37) 
with Cg an exchange coefficient and qo(z ) the humidity at ocean surface. Assuming the air is saturated 
with humidity at ocean surface, qo(z ) is estimated from sea surface temperature by: 


= 5107.4 
qo(Z ) = sat = p = 


0o(z7) 


Evaporation flux E allows the definition of latent heat flux Qr in (31) as: 








8 
640380 exp | 


oO 


| (Large 2006) (38) 


Qr=AÂE (39) 
with A = 2.5 x 108 J.kg”!. 


A.2.4 Fresh water fluxes 


Lastly, one needs to define interface conditions for salinity S,, the latter being described by: 





So ô | OSo | 
= K? 40 
Ot  OzV ) Oz em) 
The interface condition reads: 
K?ð-So = F (41) 
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where F is the fresh water flux, so that: 


with P the precipitation and E the evaporation (given by equation 37). Constant 0.001 being needed to 


get the right unit. 


F=0.001(E- P)So(z ) 


A.2.5 Summary of interface fluxes 


Momentum flux 
Sensible heat flux 
Evaporation 
Latent heat flux 
Fresh water flux 
Solar radiative flux 
radiative fluxes 
Net heat flux 


Qr = AE 


Qs, 
Qr, QLI 


A.3 Summary of continuous equations 


F=0.001(E- P)So(z ) 


The single column model equations can be summarised as: 

































































Qnet = (1 — B)Qs, + Qr, + QL) + Qe + Qn 





ôua ô | a Ua ðuo 0 | o 0Uo 
zc = = =—fk +—|K 
3: fkx (us uag) + 3- [ss az | t fkxu az V" az 
00, ð| .,00, 005 O0 |,5,000, Qs 
=: — —. = —| KT — 
at az [s dz Jes (a 915) 0t Oz dz” poC? 
Oda ô a0da So ô | 1 
= — = K 
ðt dz (x: az | ^s (da= 415) üt Oz dz 
a Ua T. o do + 
™ gz B Pa " Oz Po 
09a = QH o 000 Qs = Qnet 
* de pals ' üz  poCp  poCp 
Os 
eae = E T 2 = 0.001 (E — P) So(z ) 
a 
Ou, = OU, -0 
Oz Oz 
00, 00, -— 000 (z, 0) 
Az zi Oz 3 OZ. zs 
Oda _ 0S, = 0S, (z,0) 
Oz | dz ` Óz  |z.ho 
Atmosphere model Ocean model 


e K”? and Kj? are estimated through turbulent closure schemes TKE (Cuxart et al. 2000) and KPP 


(Large et al. 1994) for tracers and momentum respectively 


* T, Qnet and F = E — P are estimated from bulk formulae 
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| 


(42) 


T = paCa || ua (z^) ^ uo(z ) || (ua(z*) - uo (z )) 
Qu = paCy | ua(z*)) - uc(z 7) Il (8a(z*)) - 05 (27) 
E= paCz || Ug(Z*) —Uo(z ) ll (qa(z*) — Gsat) 


on Og x [0, T] 
on Og x [0, T] 
on Og x [0, T] 
on I x [0, T] 

on I x [0, T] 

on I x (0, T] 

on ang x [0, T] 
on ang" x [0, T] 


on ang" x [0, T] 


(43) 


* Ocean atmosphere coupling at the interface is done through turbulent fluxes 


Even though this is simplified model compared to the complete 3D primitive equations, it allows to 
study two major problems of the operational models, namely the computation of the interface flows via 
the bulk formulations and the evaluation of the turbulence coefficient in the surface boundary layer via 
complex parametrizations. It can be a valuable tool for both the study of coupling processes and for 
evaluation of coupled data assimilation schemes 


A.4 Discretisation 


Equations 43 are solved numerically through a second order finite difference scheme in space and a Euler 
backward scheme in time. One exception being the Coriolis term in the atmospheric moment equation 
(first equation in system 43), where a forward backward scheme is used. 


A.4.1 Vertical grid 


Vertical discretisation is similar for ocean and atmosphere with higher resolution near the interface. For 
both domain Og = [0, hg] where B = a or B = o. Both vertical grids are defined through 4 parameters: 


e hg altitude of top grid level 
e Nrg number of vertical levels 
* 0,5 stretching coefficient: the more 0,5 the more the grid will be refined close to the interface I. 


e hep height of transition between uniform (for z < hca and z > hco) and non uniform grid (for z > 
Nea z < Neo) 


from these parameters one can define the altitude/depth of cells centres zy and interfaces z% +l for 
each level k. 




















B sinho 6,5 | 
(Centres) Zk = RepO + + (hg mi ee —— with k € [1, Nr] 
sinh Osa 
. (44) 
f = hepa? hacc aad ith k € [0, Nzj] 
(Interfaces) Zk+1/2 = NcpO pir + (hp — heg) sinh&,g wi »Nrp 
where of = n and o = EX. Numerical values used in this report are summarised in table 6. 
Description Atmosphere coefficients Ocean coefficients 

Number of vertical levels Ng —51 Nro = 50 
top/bottom altitude/depth ha = 2000m o = 500m 
Transition altitude/depth hca = 200m hco = 50m 
Stretching coefficients m2 Oso = 6.5 





Table 6: vertical grids coefficients for both models 


24 


A.4.2 Time discretisation and coupling strategy 


For the sake of simplicity, both models use the same time steps (At; = Ato = At = 30s). This is probably 
not an optimal choice for a testbed, since it avoid one important difficulty in air-sea coupling. It should 
be revised in future versions of the code. 

In order to describe the coupling between ocean and atmosphere, let x; = (ua, 0a, qa)” and x, = 
(uo, 80, So)! the atmosphere and ocean states. Equations (43) can be represented by: 


MaKa) Mo (Xo) on Qg x [0, T] 

«64 (Xa) = Foa «€, (Xo) = Foa on Tl x [0, T] (45) 
OXa — OXo _ eu 
zz B Pru on 09; x [0, T] 


Implicit formulation in time implies that a simple exchange of fluxes between domain is not enough 
to find a consistent solution to equation 45. So unless the implicit problem is solved all at once for 
both model altogether, which is impracticable for realistic applications, one needs to use an iterative 
algorithme such as the Schwarz Waveform Relaxation described in section 2. Here we chose the multi- 
plicative form this algorithm: 





Ma”) Mo(x*) on Qao x [0, T] 
*€,xb-zk! Gai FE, on T x [0, T] (46) 
k k 
ox p, % B,  on0Q%* x [0, T] 
Oz Oz 


where k represents the SWR iterations, that runs until the convergence criterion is met: 


k+1 k 
oa 4515€ 


It can be transformed into a more classical asynchronous coupling by performing only one iteration. 


A.5 Numerical implementation 


The coupled SCM is made of 70 000 lines of Fortran/C++ including tangent and adjoint and interfaces. 
SCM tangent and adjoint models were generated using the automatic differentiation tool TAPENADE?’. 
It has been interfaced with OOPS for data assimilation use and will then be made available to the com- 
munity and be proposed as a reference test case to be included in OOPS. For now, it is considered as a 
single model by OOPS, which is suboptimal and will be revised. However it has allowed to implement all 
assimilation algorithm presented in this report. 


B Single column model reference test case 
This appendix describes experimental settings that are used in section 6. The single column model pre- 


sented in appendix A is set up so that it mimics the behaviour of a mid-latitude air-sea column (coriolis 
parameter f = 1074 s71). 





?http://www-sop.inria.fr/tropics/tapenade.html 
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B.1 Atmosphereinitial conditions 

For the atmosphere, initial conditions are 
Ug(Z, t = 0) = Ugo = 15m-s! 

1 


ZEOg 
Valz, t = 0) = Vao =3 m-s 


NZeret (47) 
— — zK  zeQ, withN=0.01s ! and 6 = 286K 


zEeQ, 
6a(z, t = 0) = 049 = OT + 


qa(z, t = 0) = qao = 0.01 kg/kg ZEQa 


where N is the Brunt-Väisälä frequency and g is the acceleration of gravity constant. 
While for the ocean column, one sets up 


1 


Uo(z,t=0)=0m:s zeQ, 
Vo(z,t=0)=0m-s! zeQ, 
=0= .10 3. : Z ag. A (48) 
6,(z, t=0) 2 5--1.475-10 3. 24.5 eph) +8 exp |z) K zeQ 
So(z,t =0) = S=35- (1.75: 107*- 20.8: exp(-—) +0.2-exp(—}) psu ZEOg 
200 400 


The system then undergoes a 10-day spin-up, which leads to reference initial vertical profiles showed 
in figure (7) for both atmosphere (top) and ocean (bottom) variables. 


B.2 Large scale relaxation 


The atmosphere model includes some large scale relaxation terms uag, O75 and qrs (see equations 43). 
They are defined so that the system shows enough variability, in particular so that it switches from stable 
to unstable regimes on a short enough time window. In addition, one would like to observe a diurnal 
cycle for the temperature. 

This is achieved setting: 








Ugo + Uup: l(z, t) siz < hii 
UaG(Z, t) = is d . m zeQOQg, te [0, T] 
Uao SiZ > him 
Vao + Vun: l(z, t) siz < hi 
vaclz, =] ^" ts z€ Qs, t € [0, T] 
Vao si z > him (49) 
5 27t Qs(0, 7) Pd 
Ors(z, t) -0 COS . +1.l-ex zeQ,,tel0o,T 
LL G | = max (Qs(0, 0) Fe : 
qis(Z, t) = dao ZEQzg, te [0, T] 





where T, = 86400 s, U, = 15 m.s !, V, = 13m.s ! and: 


NZ 
l(z, o-ew[-5 5 }-sin( 7 | 
i lim 








with u = 25-104 s, jim = 400 m and ø = 6- 10^ s. Function / along with Uup and V, allows to define 
a wind variation for 0 < z < hjim, which looks like a Gaussian, centred around the 3' day of the time 
window. Large scale relaxation term 0, is defined so that it mimic the diurnal cycle. Figure (8) shows 
these large scale forcings on a 7 day window. 
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Figure 7: Initial vertical profiles for atmosphere (top) and ocean (bottom) reference test case. (a) shows 
wind u,(z,0) and v,(z,0) profiles and (b) atmosphere temperature and humidity 0,(z,0) and humidity 
qa(z,0) profiles. (a) shows ocean current velocity profiles u5(z,0) and v,(z,0) for the first 50m and (b) 
shows ocean temperature 0,(z,0) and salinity So(z, 0) profiles. 


B.3 Prescribed ocean-atmosphere Fluxes 


As mentioned before, Qs,, Qr, and P fluxes are exogenous of our system and have to be prescribed 


To do so, one first needs to define the surface pressure 


p(0, t) = 101320 Pa Vt €[0, T] (50) 


Solar flux Qs(0, t) is then defined so that both a diurnal cycle and a day-to-day variability are present 
(see figure 9a). 
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Figure 8: Large scale atmospheric forcings. (a) and (b) show the zonal and meridional geostrophic winds 
and (c) shows the large scale potential temperature. 





ee) Vt € [21600 + k- /,64800 + j- J] 
20? with J = 86400s and j € [0,6] 
0 otherwise 


Var(£) : Qrnax * EXP |- 
So — 





ith Var(t) = 0.8 + 0.2- 
wi ar(t) cos [cca] (51) 


Qmax = 420 W.m ? 
o = 7200s 
u=43200+j-Js 
There is then no solar flux at night (between 18 h and 6 h) and it follows a Gaussian-type evolution be- 


tween 6h and 18h. Daily maximum is set via Var(t) which adds a day-to-day variability throughout the 
time window. 
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Precipitation flux P(t) is set so that strong precipitations occurs at the same time as the strong geostrophic 
wind event mentioned in the previous section (around 3" day, see figure 9b). This can be seen as the 
signature of a depression passing through our atmosphere column. 





EN 
PO = ca Sew [— p | mts? Vt € [0, T] 


0.001 202 m 
with u= 200000 s me 
a =30000s 








P(t) (cm.day-!) 


Qs (t) (W.m7?) 




















1 2 3 4 5 6 7 ^ z = 
t (day) t (day) 





(a) Prescribed solar flux Qs(0, t) (b) Prescribed precipitation P(t) 
Figure 9: Reference prescribed solar flux at ocean surface Qs(0, t) (a) and precipitation P(t) (b) over time 


Lastly waveband radiation flux Qr, is set constant throughout the time window. 


Qr (t) 2350W.m ^ Vtel0,T] (53) 


B.4 Modelsimulations 


Figures (10) show profiles of wind and current velocities over all the time window. We can notice an 
increase responding to that of geostrophic winds. This increase had the effect of raising the level of the 
atmospheric boundary layer around 1000 m before returning to its previous altitude. We can also observe 
that the air temperature varies according to the diurnal cycle. This comes from changes in ocean surface 
temperature directly impacted by solar flux Qs. Figure (12) represents temperature differences at the air- 
sea interface between both media. This diagnosis allows to determine if the system is in a stable (0, > 00) 
or unstable (8, < 0o) regime. 
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Figure 10: Colors represent zonal (a,c) and meridional (b,d) atmosphere wind and ocean current veloci- 
ties components and black isolines represent temperatures. 
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Figure 11: First guess for ocean interface conditions for SSU and SSV (a) and SST b 
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Figure 12: Differences between atmosphere and ocean surface temperature 04(z*, t) ? 0o(z*, t) 


31 


C Elements of convergence study 


Let us recall that the aim of any data assimilation methods mentioned in this report is to minimise, albeit 
approximatively, the following cost function: 


1 T 1 
Jo) = 5 (x-x^) p! (x-x!] +5 (AM) -y°) R^! (Z6 C490) - y?) (54) 


where M is the coupled model. 

Since iterative schemes are used for solving the optimisation problem it is important to study their 
convergence. To that matter, we use a different classification than that of the core of the report. On the 
one hand we consider incremental and approximate incremental 4DVar schemes and on the other hand 
splitting-based 4DVar. Following this classification for the algorithm presented in this report, FCM can 
be seen as a direct application of incremental 4D-Var, while PCM is an approximate incremental 4DVar. 
WCM is a splitting based 4DVar and CERA can either be seen as a splitting 4DVar or an approximate 
incremental 4DVar. 


C.1 Incremental and approximate incremental 4DVar 


Incremental 4DVar is a well known and thoroughly studied variational data assimilation scheme. In 
particular Gratton et al. (2007) (GLN2007 hereafter) studies convergence properties for approximate in- 
cremental 4D-Var, where the tangent linear model used in the inner loop is an approximation of that 
of the original problem. In this section we shamelessly exploit results from GLN2007 that are directly 
applicable to our problems. 

First let us reformulate the 4Dvar problem in a form compatible with GLN2007 and recall their main 
results. If n is the size of the control vector and p the size of the observation vector, let F : R” — R”+P 
such that 


B2 (x -x^) 
F(x) =| R-!2( 7M) y^) (55) 
Equation 54 can then be rewritten 
l 2 
Jo) = 2 I.FG91l5 (56) 
B2 
Denoting Fx = | RBM | the jacobian (tangent linear) of F differentiated around x, gradient 
XX 


and Hessian of J read 


VxJ=FĪF&) e m" (57) 
VVxJ-FIF,-Q() € m"*" (58) 


where Q(x) denotes the second order terms 

Incremental 4DVar, which can be seen as a Gauss-Newton algorithm, solves problem 55 by successive 
approximations neglecting the second order terms. 

Theorem 4 of GLN2007 states that a sufficient condition for this algorithm to converge to the mini- 
mum of the original problem is that there exist a sequence n < 1 such that 








=] 
| Qa) (FoF) | <n (59) 
2 
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Algorithm 1: Incremental 4D-Var 


while not converged do 
compute F(x“) and differentiate; 
solve inner problem Le Fu dx = -F7 F (x); 
update reference x **P =x + 5x (0 

end 


Meaning that the change of curvature of J must not be too strong. Additionally, for convergence to be 
ensured the first guess of the algorithm x^ has to be 'close enough' to the global optimum (depending on 
the aforementioned curvature). 

In general, due to non differentiabilities in the model and in order to save computing time, the ja- 
cobian used in the inner loop is only an approximation of F,w. For instance PCM, with only part of the 
coupling being accounted for, is a good example for such an approximation. CERA goes a little bit fur- 
ther, neglecting the coupling processes in the inner loop altogether. In order to describe this approximate 
incremental 4DVar (a.k.a Perturbed Gauss-Newton in GLN2007), let us denote F, this approximate tan- 
gent operator. 


Algorithm 2: Approx-Incremental 4D-Var 
while not converged do 
compute F (x) and differentiate; 
solve inner problem ES E, dx = -FF (x); 


update reference x/** =x + ôx 


end 


Then GLN2007’s theorem 6 tells us that the sufficient condition becomes 


I = (Be ES Qo) (Er Ee) 








<n® (60) 
2 


Note that if F,« =F, this is equivalent to condition 59. Condition 60 is a bit more difficult to interpret 
in a geometrical way. However it states that for convergence to be guaranteed, the requirement is that 
FIF, has to be a good approximation of FIF, + Q(x) (the Hessian of the approximate inner problem). 

Even if the Approximate 4DVar converges, there is no guarantee that it will converge toward the same 
minimum as the original problem. In general it is not the case, and theorem 7 from GLN2007 gives a 
superior bound of this error at optimum 


Iz -x' Ios — | (FE - Ej) rl = — [e Fl (6) 
1-v 1-v 
where F* = (FTP) F7 denotes the Moore Penrose inverse and 0 < v < 1 is the upper bound of the con- 
vergence speed of the non approximated problem and depends on second order terms Q (v = 0 in the 
linear case). In vernacular words, this can be interpreted as the less regular / is, the less the inner model 
can be approximated. 
At the fixed point, the perturbed Jacobian F,« must, be such that FI F(x*) 
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Remark: This section only focuses on convergence (or not) ofthe outer problem. The inner problem is 
quadratic, so it always converges. Regarding the speed of convergence, for outer problems it is directly 
driven by 74 values (the smaller the better), while for inner problems it depends on the conditioning of 
F7 th) F,w. However CERA can perform the minimisation separately on each media, so the inner problem 
is ‘likely to converge faster, depending on the conditioning of the sub problems. 


C.2 Splitting methods 


WCM, and CERA to a certain extent, can be seen as a splitting algorithm. Indeed, the original problem is 
split into two somewhat independent problems. In particular WCM can be seen as resembling a parallel 
version of the antique Alternating Direction Method of Multiplier (Fortin & Glowinski 1985, ADMM), 
which is used in machine learning. 

ADMM is meant for linear problems where the minimisation problem can be rewritten into 2 inde- 
pendent constrained minimisation problems e.g. finding 


argmin J (x) 


is equivalent to finding 


argmin (Jg(Xq) + Jo (Ko) 
under constraint Ax; = Bx, 


The ADMM algorithm is then: 


Algorithm 3: ADMM 


initialise x = x? ,x® = x? and c® =0; 


while not converged do 
x^? = argmin,, (Jala) + yllAxa + Bx +e |); 


xD argmin,, (o9) + ylAx ^P + Bx, + c I"); 


c ^D = eO + Ax 5D - px +0 
end 


ADMM is known to converge if J, Ja and J, are proper, convex and lower semi-continuous functions, 
which is generally the case in data assimilation if 7 and M are linear. 

In CERA, which is not really an ADMM, y = 0 but J, and J, are updated at each iteration. WCM on 
the other hand is quite close, y||Ax; + Bx; + cl? resembling the /? term. 

In the non linear case, in order to ensure convergence, one could imagine a 3-level algorithm where 
the inner problem of a non approximated incremental 4DVar would be solved through an ADMM al- 
gorithm. Such an algorithm would require the same conditions as for incremental 4DVar to converge. 
However, it would be quite expensive since it would multiply the number of minimisations. WCM pro- 
poses to squeeze it into a two-level minimisation, just like incremental 4DVar, where both non linearity 
and /? update are done at the same time (during the outer iteration). Convergence study of this variant 
remains to be done. 
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C.3 Discussion 


Considering coupled data assimilation schemes as (perturbed) Gauss-Newton algorithms allows for reusing 
the theoretical convergence results. However in the non linear case, theory only provide sufficient condi- 
tions (i.e. the algorithm may converge even if they are violated). Moreover these conditions are difficult 
to evaluate due to the second order terms. On the other hand, in the linear case these sufficient condi- 
tions become also necessary and are far easier to evaluate if tangent and adjoint models are available. 
Consequently, if full convergence of the algorithm were the main concern, the best option would proba- 
bly be to use three level strategy where the inner problem would itself be solved through either Perturbed 
Gauss-Newton or splitting. In practice full convergence is seldom sought for, and the number of iteration 
being limited, speed of (partial) convergence at the early stage of the algorithm is the most important. 
Theory also gives a measure of inner convergence speed (the nxs), but yet again, this can only be really 
exploited in the linear case. Evidence actually shows that ECMWF operational uncoupled incremental 
4D-Var does not converge as outer loops is concerned but it is very efficient at providing a good analysis 
with a limited number of outer iterations (Trémolet 2007). 
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